An efficient finite element method appfied to quantum biUiard systems 
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An efficient finite element method (FEM) for calculating eigenvalues and eigenfunctions of quan- 
tum billiard systems is presented. We consider the FEM based on triangular Ci continuity quartic 
interpolation. Various shapes of quantum billiards including an integrable unit circle are treated. 
The numerical results show that the applied method provides accurate set of eigenvalues exceeding 
a thousand levels for any shape of quantum billiards on a personal computer. Comparison with 
the results from the FEM based on well-known Co continuity quadratic interpolation proves the 
efficiency of the method. 
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I. INTRODUCTION 

There has been much interest in characterizing the 
quantum manifestation of classically chaotic systems 
[U, [3] , since McDonald and Kaufman's pioneering investi- 
gation on the statistical characteristics of eigenvalues and 
eigenfunctions [3j]. The quantum billiard, which is rep- 
resented by the two dimensional stationary Schrodinger 
equation of free particle with satisfying the Dirichlet 
boundary condition, is an intensively studied model sys- 
tem in the field of quantum chaos due to its simplic- 
ity. The integrability of corresponding classical billiard 
depends solely upon the geometry of boundary. The 
quantum billiard can be also expressed by the scalar 
Helmholtz equation, for example, which describes the 
electromagnetic field inside a flat microwave resonator. 
In that context, microwave experiments played the role 
of analog computation of eigenstates in quantum billiards 

There are several numerical methods, which have been 
dominantly adopted by colleagues in this field, for cal- 
culating eigenvalues and eigenfunctions of the quantum 
billiards such as the boundary integral method (BIM) 



(reviewed in Ref. 
method (PWDM) 



the plane wave decomposition 



cornposr 
the scaling method [1, , and 
the conformal mapping method [l^.Tl3l. . Let us briefly 
review the mentioned methods. The BIM, which is a rig- 
orously established method, reduces the problem of two 
dimensional stationary Schrodinger equation to a one di- 
mensional integral equation. In result, each root of the 
Fredholm determinant constitutes the set of eigenvalues. 
In practice, the determinant does not become zero due 
to the discretization error and the BIM approximates the 
minima of the lowest singular values to the eigenvalues. 
Though many successful applications, the BIM has some 
shortcomings. One of them is that the calculated results 
can include additional, i.e., spurious solutions which cor- 
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respond to roots of outside scattering problem with Neu- 
mann boundary condition (isj . For non-convex geometric 
billiard, the detection of spurious solutions is not a simple 
numerical task and the BIM failed in the isospectral 
drum introduced by Gordon, Webb, and Wolpert due 
to its strong non-convexity [13] . The feasibility of miss- 
ing eigenvalues is an another weakness of the BIM. For 
higher lying eigenvalues, the spectra become more denser 
and the correct detection of minima is a serious problem. 

The PWDM, which has been introduced by Heller 0], 
is a rather heuristic method in the context of quantum 
chaos. It is appropriate for computing high lying eigen- 
states but incongruent for studying of spectral statistics, 
because only a few selected eigenstates can be calculated 
with many intermediate missing. Also the PWDM can 
fail in non-convex or multiply connected (e.g., containing 
a hole) billiards [3]. In the literatures, there has been 
a considerably efficient numerical method, that is, the 
scaling method derived by Vergini and Saraceno [9|. It 
represents the boundary norm as a function of energy by 
the use of scaling. In result, the authors of Refs. [lO] 
obtained all eigenvalues (without any missing) within a 
narrow energy range, which lie close to a chosen reference 
value, in a single computational step. The efficiency of 
scaling method is obvious from that the BIM can locate a 
single eigenvalue in a single computational step. For spe- 
cific geometric billiard for which the conformal mapping 
onto the unit disk is sufficiently simple (e.g., so-called 
Robnik billiard \^M) , the conformal mapping method de- 
rived by Robnik [l2| have provided accurate set of eigen- 
values jlTj l . Recently, new approaches that combine each 
ideas of above mentioned methods have been studied, for 
the BIM and the PWDM ^ and for the BIM and the 
scaling method . Concerning the scattering quantiza- 
tion method, an efficient improvement has been carried 
out in Ref. [2l|. 

The finite element method (FEM) is one of the most 
widely accepted numerical methods for partial differen- 
tial equations in various fields of science and engineering 
[12, [H, [131 . Compared with previously mentioned meth- 
ods, the FEM has obvious advantages that it has almost 
no limitation on the geometric complexity of billiard (see 
the results in Ref. (25|) and provides in a single compu- 
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tational step a set of all eigenvalues and eigenfunctions 
up to maximal level allowed by memory allocation. How- 
ever, the accurate computation of high-lying eigenstates 
using the FEM is conventionally more difficult than the 
case of other mentioned methods, since the FEM dis- 
cretizes not only the boundary but the whole domain 
of billiard (it needs more memory storage). Thereby, 
though its obvious advantages, the FEM has been ap- 
parently overlooked in the field of quantum chaos. As 
far as we know, there have been only a few studies 
[lllil,!!^,!!!!!!] where the FEM is used for calculating 
the eigenstates of quantum billiards. Among those stud- 
ies, Heuveline showed an effective FEM only requiring 
0(N) memory allocation by usin g th e p-finite elements 
basis and the sparsity of matrices [28|, l29| . Note that the 
FEM commonly gives rise to sparse matrices but usual 
FEMs do not take advantage of the sparsity (it contains 
quite difficult numerical tasks) and need 0{N'^) memory 
storage, where N is the number of total nodes. 

The aim of this paper is to show the validity of FEM 
as a numerical method for calculating eigenvalues and 
eigenfunctions of quantum billiards. For that purpose, 
we present an efficient FEM based on the Hermite in- 
terpolation. In each element, the wave function is in- 
terpolated by quartic polynomials involving nodal val- 
ues of wave function and its first derivatives, namely, the 
adopted interpolation basis admits the Ci continuity. By 
applying the method, we calculate the eigenvalue spec- 
tra of unit disk (integrable), the Robnik billiard (convex 
geometric chaotic), and the spiral-shaped billiard (non- 
convex geometric chaotic). We show that the method 
provides accurate set of eigenvalues exceeding a thousand 
levels for any shape of quantum billiards on a personal 
computer. Comparison with the results from the FEM 
based on well-known Co continuity quadratic interpola- 
tion proves the efficiency of the applied method. Note 
that, by virtue of the Ci continuity, the method handles 
well problem that treats values of first derivatives of wave 
function at the boundary such as the Neumann boundary 
condition. 

The rest of paper is organized as follows. In Sec. II, 
we outline numerical procedures of the FEM based on Ci 
continuity quartic interpolation. The results of numeri- 
cally calculated eigenvalues and the analysis of spectral 
statistics for unit disk, the Robnik billiard, and the spiral- 
shaped billiard are presented in Sec. III. Conclusions are 
given in Sec. IV. 



natural units h = 2m = 1. 

The first step of applying FEM is that the domain of 
billiard is discretized into finite elements, i.e., mesh 
generation. The shape of element and the number of 
nodes in each element are determined according to the 
type of interpolation polynomials (i.e., shape functions). 
Here we consider the FEM based on the triangular Ci 
continuity quartic interpolation, which has been derived 
by Specht [30| and known that passes all patch test, i.e., 
a condition for assessing FEM convergence for arbitrary 
mesh configurations (see Chapter 11 of Ref. [131 )• 
cordingly the domain of billiard is discretized into trian- 
gular elements. In each element, there exist three nodes 
that locate at vertices of triangle and each node has three 
degrees of freedom correspond to wave function and its 
first derivatives {ip, dtp/dx, and dip/dy). Thereby each 
element has actually nine nodes. In each element the un- 
known function, i.e., wave function ')jj{x, y) is represented 
as a linear combination of shape functions multiplied by 
as-yet- unknown nodal values of wave function and its first 
derivatives. The shape function is defined only over a 
given element and has zero value at outside of it. An ex- 
plicit representation of shape functions will be postponed 
for a while. The numerical procedure of FEM requires 
that each node has three indices; local index il, element 
index ie, and global index i. The mesh generation com- 
pletes the mapping from local and element index {ilth 
node of ieth element) to global index {ith global node). 

On numerical calculations in this paper, we use two 
considerable and freely available mesh generators; the 
DistMesh [IH and the Triangle The DistMesh is a 
Matlab based mesh generator that finds node locations 
settled down a equilibrium state in a truss structure. A 
geometry of domain is represented by the signed distance 
function from node to closest boundary dU,, negative in- 
side the domain. It generates high quality meshes, i.e., 
almost equilateral triangles but can be faced with diffi- 
culty for complex geometric boundary. The Triangle is 
a robust Delaunay refinement code. The user-supplied 
data, which contains the information of nodes placing on 
the boundary, are employed for specifying the domain of 
billiard. The quality of meshes is controlled by the con- 
straint of minimum angle (up to 34 degree) and maximum 
size of triangle. The Triangle has almost no limitation on 
the complexity of geometry. 

Now we can take the next step. Eq. (1) can be ob- 
tained from the condition that the action 



II. NUMERICAL PROCEDURE 

The quantum billiard is governed by two-dimensional 
stationary Schrodinger equation of free particle 

- V^VW = E^ir), for r e VL (1) 

with satisfying the Dirichlet boundary condition ip{r) = 
at the boundary of domain Note that we use the 



j ds{\/ip*[r) ■V^{r)- E^*{r)7P{r)) (2) 

is minimized with respect to variation of 'tp*{r). Here 
■!/)(r) and tp*(f) ^re considered to be two independent 
variables. Then, the action integral of Eq. (2) is dis- 
cretized into integrations over each element as 

ne 
ie 
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FIG. 1: A triangular real element and its mapping onto the 
parent element, i.e, a right isosceles triangle. 



where ne is the number of total elements. In Cartesian 
coordinates, the discretized action integral A^'^^ is repre- 
sented by 



-E 



dx dx dy dy' 



dxdy 



(4) 



P = { 1 - C - r?, (1 - ^ - ri% ^v, - ^ - V), 

(3(1 - /i3)(l - ^ - ry) - (1 + 3/X3)C + (1 + 3^3)7?) , 

ev + ^il-C-rMv (8) 
(3(1 - - (1 + 3Mi)r/ + (1 + 3/xi)(l - ^ - , 

(3(1 - H2)V - (1 + 3m2)(1 - ^ - r?) + (1 + 3^2)0 } 

where /i,,, = (Zc - lb)/la and ?a = (a;6 - Xc)^ + {yb - ycf' 
for the cyclic permutation of a, h, c. The coordinate 
transformation between real and parent elements is given 
by the following Co continuity linear shape functions 



ar = ^xf)iV,(e,,7), 2/ = ^2/f'iV,(C,r,). (9) 



where ijj = tp{x,y) and tp* = tp*{x,y). For simple com- 
putation of A'^^'^^ , it is advisable that the integral domain 
of each clement is transformed into a regularized domain 
(called parent element) as 



= dxdyf{x,y)= d^drj g{^,T]). (5) 

J J ieth J J parent 



Figure 1 shows a transformation of triangular real ele- 
ment into the parent element. The quantities in Eq. (4) 
are altered into (C,*?) notations as the foUowings. The 
wave function ipi^^v) is interpolated by the Ci continu- 
ity quartic shape functions {Hii{^,ri)} multiplied by as- 
yet-unknown nodal values of wave function and its first 
derivatives {tpu} as 



where Ni{^, r?) = l-^-rj. A"2(^, q) = ^, and N3{^, if) = r). 
Then, the infinitesimal surface dxdy is represented by 



''^^^(ll-ll)*^'- « 

The partial derivatives of wave function are transformed 
into (^, rf) notations as follows: 



dip{x,y) d 



dx 



dx 



d(, dx 



+ 



drj dx 



(11) 



(6) 



d^{x,y) _ d I s^jUe) 



dy 



dy 



dx dy dx dy dx dy 

For simple representation, we use a cyclic property of 

the applied shape functions and introduce the following 
notation Ha = {i?3a-2(^, -ff3a-i(^, J?), -f^3a(^, J?)} for 
a = 1, 2, 3. Then, the shape functions are represented by 



H 



- Pa+3 + Pc+3 + 2(Pa+6 " Pc+e) 
mb{Pc+6 - Pc+s) + mcPa+6 

-nb{Pc+6 - Pc+3) - ncPa+e 



(7) 



where nia = Xc—Xb, ria = yb—Vc-, and a, 6, c are the cyclic 
permutations of 1, 2, 3. The nine polynomials {Pi{^,rf)} 
in Eq. (7) are expressed as 



a ^ 



7,(ze) ( dHii{^,ri) dS, _^ dHu{i,ri) drj 



d$, dy 



df] dy 



By applying Eqs. (6)-(ll), the discretized action inte- 
gral A^^^^ is represented by 



il,jl 







1 .1-7, 

dr] I d^ /(C, ri)ilS 1 



US 







1-77 



(ie) 



dn j d^g{^,r^)us]'4^f\^'^) 



where /(^, 77) and g(^, rj) arc sixth and eighth order poly- 
nomials, respectively. By using the optimized quadrature 
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rule over the triangle, which has been derived by Duna- 
vant [33|, Eq. (12) can be exactly integrated as 

dv f " /(e, ^) = E ■ (13) 

Jo Jo j 

where (^j, ry^) is a quadrature point and uji is a weight. M 
equals 12 and 16 for sixth and eighth order polynomials, 
respectively (see the table in Ref. (s^). Then we obtain 
the following result 

iljl 

Now we add up the discretized action integral of Eq. 
(14) according to Eq. (3). Then the action is represented 
as 



N 

A = Ev;*.(L,,,-ii;M,,,)-v;, (15) 

where Li_j (Mi^) is a summation of all L^l'^j^ fo^' 
which satisfies that ilt^^ and jl^^ nodes of icth element are 
mapped into ith and jth global nodes, respectively. N is 
the number of total global nodes. In our case, it is equal 
to three times the number of physical global nodes, since 
each node has three degrees of freedom. 

The adaptation of Dirichlet boundary condition is im- 
plemented as follows. If the kth global node places at the 
boundary of billiard dVL^ the nodal value of wave function 
ipk equals zero. It requires that the entries of kth column 
of L and M matrices become zero and also the entries 
of kth row set to zero since ip1 = Q. In practice, this is 
achieved by dropping the kth row and column. Then, the 
dimension of L and M matrices are reduced by — A^f, 
where Ni, is the number of global nodes located at d^. 

Now we vary the action A with respect to nodal values 
"0* and invoke the principle of least action. In result, we 
obtain the discretized version of stationary Schrodinger 
equation in the form of generalized eigenvalue problem 



L* = £;A/* (16) 

where L and M are N — Nf, dimensional real and sym- 
metric matrices and is a column matrix, which consists 
of nodal values Tp. We use the LAPACK routine DSPGV 
[3^ for solving the generalized eigenvalue problem. Note 
that this routine requires about 2{N — N^f memory stor- 
age. Finally we obtain the set oi N — Ni, eigenvalues {Ei} 

and eigenfunctions {vl/i} where "^J = ^'tpi ■ • • "ipN^ 
through restoring zero nodal values at nodes on dfl. 



III. NUMERICAL RESULTS 
A. The unit circle billiard 

We firstly consider an integrable billiard whose bound- 
ary is given by the unit circle for testing an efficiency of 
the FEM presented in Sec. II. In this case, the eigenvalue 
spectra are exactly known and given by the sorted set of 
/i^j. where iJ,jk is the kth root of jth Bessel function of the 
first kind with considering the degenerate case on j 7^ 
for j = 0, 1,2, • • • , and A; € N. 

We investigate the relative error between exact and 
numerically calculated eigenvalues 



£rel{i) = ^^V^ (17) 

where {Ei} is the set of exact eigenvalues. For the results 
obtained from the FEM based on well-known tri ang ular 
Co continuity quadratic interpolation (see Refs. [U, H^l 
for its shape function and isoparametric transformation), 
we also compute the relative error Srei- In Fig. 2(a), we 
plot the relative error Srei for both cases of interpolation 
basis. For the C'l continuity quartic interpolation, we 
calculate 21500 (equals to TV — Ni,) eigenvalues and the 
result of relative error is drawn by a black line. For the Co 
continuity quadratic interpolation, we obtain 21647 (al- 
most same number of the above) eigenvalues and depict 
the relative error as a orange line. Note that all calcula- 
tions in this paper are performed on a personal computer 
possessing 2.4 GHz Quad-Core CPU and 8 GByte mem- 
ory. So the available maximal number of level allowed by 
memory allocation is limited about 22000. In both inter- 
polations, the relative error Srei increases as the number 
of level n increases. However, the method based on the 
Co continuity quadratic interpolation seriously loses its 
accuracy after a few hundreds levels. The results of Fig. 
2(a) proves an efficiency of the FEM based on the Ci 
continuity quartic interpolation. 

We consider the numerical result for which the relative 
error Srei is smaller than 3 x 10"'^ as accurate eigenvalue. 
In Fig. 2(b), the relative error for the Ci continuity 
quartic basis is depicted up to 1500 level and the result 
is smaller than 3 x 10^'^ in this range. Then we obtain 
accurate set of 1500 eigenvalues for the unit circle billiard 
with this criterion. The foUowings will show that above 
conjecture is reasonable. 

Another method for testing the accuracy of obtained 
results is checking out whether the eigenvalue spectra 
are complete without any intermediate missing. For the 
system where the analytic eigenvalues are not available, 
such method has no alternative. It can be performed by 
investigating the spectral staircase function N{E), which 
counts the number of energy levels below E. The spectral 
staircase function can be divided into a smooth and a 
fluctuating part 
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FIG. 2: (Color online) (a) The relative error Srei for the Ci continuity quartic interpolation (the Co continuity quadratic 
interpolation) is drawn by a black (orange) line, (b) The relative error £rei up to 1500 level for the Ci continuity quartic basis, 
(c) The 5„ up to 1500 level for the Ci continuity quartic basis. 



(18) 



The smooth part N{E) is represented by the gcnerahzed 
Weyl's law ^ 



N{E) 



A L r- 

47r 47r 



C 



(19) 



where minus and plus sign correspond to the Dirichlet 
and the Neumann boundary condition, respectively. A is 
a area of billiard, L is a length of the perimeter, and C is 
a correction constant for the curvature and corners given 

by 



C 



1 

12^ 



k{s) ds 



24 ^ 



TT 



(20) 



with local curvature k{s) and ith corner angle a.;. 

The so-called 5n quantity, which is equivalent to 
Nfiuc{En), is a good measure for the completeness of 
obtained results 



N{E,,) 



(21) 



where N{En) := n — i. For complete eigenvalue spectra, 
it has been well known that (5„ fluctuates around zero. 
In Fig. 2(c), 5n is drawn up to 1500 level for the Ci 
continuity quartic basis. It certainly fluctuates around 
zero and shows that the obtained eigenvalues are com- 
plete up to 1500 level. Figure 2(c) also proves that the 
above criterion for the relative error is acceptable. 



The Robnik billiard 



billiard is defined by a quadratic conformal mapping from 
the unit circle 



^ X + iy = z + \z 



(22) 



where z = e**", (j) G [0, 27r], and A e [0, i]. With increas- 
ing A from zero, the boundary is continuously deformed 
from the unit circle. At A = i, the mapping of Eq. (22) 
is no longer conformal and the billiard has a cusp. Such 
limit case of the Robnik billiard is also called the car- 
dioid billiard. It has a symmetry line at y = and the 
desymmetrized billiard is twofold; odd and even symme- 
try satisfying the Dirichlet and the Neumann boundary 
condition at the symmetry line, respectively. In Fig. 3, 
we show the cardioid billiard and its desymmetrized ver- 
sion. 

It has been proven that the cardioid billiard is ergodic, 
mixing, and a K system, i.e., fully chaotic system [36| and 
its spectral statistics have been minutely studied in Ref. 
[stI ] . We regard the odd type of desymmetrized cardioid 
billiard as a model of convex geometric chaotic system 
and test the FEM based on the Ci continuity quartic in- 
terpolation. By applying the method, we calculate 21901 
eigenvalues of the odd symmetric case. Among those we 
obtain 1500 accurate eigenvalue spectra. As shown in 
Fig. 4(a), (5„ fluctuates around zero in this range. 



(b) 





In this section, we consider the billiard which has been 
introduced by Robnik ^ll'| . The boundary of the Robnik 



FIG. 3; The cardioid billiard, i.e., the limit case of the Robnik 
billiard; (a) full and (b) desymmetrized version. 
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FIG. 4: (Color online) For the odd type of desymmetrized cardioid billiard; (a) The (5„ up to 1500 level, (b) The cumulative 
level spacing distributions of the GOE, GUE, and the considered billiard are drawn by a black full, black dotted, and orange 
full line, respectively, (c) The spectral rigidities of the GOE, GUE, and the considered billiard are drawn by a black full, black 
dotted, and orange full line, respectively. 



We investigate two spectral statistics, that is, the 
nearest-neighbor level spacing distribution and the spec- 
tral rigidity for 1500 accurate eigenvalues. The nearest- 
neighbor level spacing distribution P(s) ds is the proba- 
bility of finding a consecutive pair of eigenstates for which 
the difference between their eigenvalues lies in the inter- 
val [s, s -l- ds]. It measures the short range correlation of 
the eigenvalue spectra. Instead of P(s), we consider the 
cumulative level spacing distribution 

I{s) = f P{s')ds' (23) 

to keep out of the binning problem about P{s). The spec- 
tral rigidity A3(L) is the mean square deviation of the 
spectral staircase function from the best fitting straight 
line over a length L, namely 

A3(L) = (min-/ de {N{E + e)-a-beY) .(24) 

\ (a, 6) L 7_L/2 / E 

It was firstly introduced by Dyson and Mehta ^] to 
describe statistics of the energy levels of many particle 
systems such as nuclei. It measures the long range corre- 
lation of the eigenvalue spectra. Through studies about 
two mentioned spectral statistics, we rescale the eigen- 
value spectra {Ei} into {E^} where E^ = N{Ei) and we 
omit the prime. After the rescaling, the eigenvalue spec- 
tra have a mean level spacing of unity and each billiard's 
own characteristic is contained on the fluctuating part 
Nfiuc{E). 

It has been widely accepted [39l| that the spectral 
statistics of classically fully chaotic systems can be well 
described by the universal laws of random matrix the- 
ory (RMT) [13]. From the RMT prediction, the spectral 
statistics are given by the distribution of the Gaussian or- 
thogonal ensemble (GOE) and the Gaussian unitary en- 
semble (GUE) for systems with and without time reversal 



symmetry, respectively. Note that the time reversal in- 
variant systems possessing specific geometric properties 
can show the GUE-like statistics. Concerned discussions 
will be addressed in next section. From Berry's semi- 
classical analysis for spectral rigidity [4l|, it has been 
also known that the universality region where the spec- 
tral statistics follow the universal RMT prediction is fi- 
nite. On paraphrasing, for fully chaotic systems, the 
spectral rigidity A3(L) shows a universal logarithmic in- 
crease following the prediction of RMT in the interval 
1 < i < Lmax- For the case of GOE, the coefficient 
of logarithm is twice that of GUE. Then, in the range 
L > Lfnax, ^3{L) reaches a non- universal saturation 
value determined by short periodic orbits of correspond- 
ing classical billiard. Lmax is called the outer energy scale 
and depends on the period of shortest periodic orbit and 
the mean level density. 

Since the cardioid billiard has the property of time 
reversal invariance, it is expected that the cumulative 
level spacing distribution I{s) follows that of GOE and 
the spectral rigidity A3 (L) is well described by the GOE 
prediction within the universality regime. In Figs. 4(b) 
and 4(c), we show the results of I{s) and A^lL) for the 
odd type of desymmetrized cardioid billiard and compare 
with those of the GOE and GUE (for numerical calcu- 
lation of I(s) and A3(L) for the RMT predictions, see 
the Ref. [33 )• As expected, the results show that the 
spectral statistics are in good agreement with the GOE 
predictions and the spectral rigidity saturates beyond the 
universality regime, which is restricted to small correla- 
tion length L. 

Note that we also calculate eigenvalue spectra for the 
even type of desymmetrized cardioid billiard by applying 
the Neumann boundary condition at the symmetry line. 
It can be easily achieved by the Ci continuity property 
of the applied shape function. We obtain the same 1500 
accurate eigenvalues as the odd symmetric case. The re- 
sults of spectral statistics for the even symmetric case are 
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well described by the GOE expectation as qualitatively 
equivalent to the odd symmetric case of Figs. 4(b) and 
4(c). We would not present these results in figure. 

C. The spiral-shaped biUiard 

In this section we consider the spiral-shaped billiard 
whose boundary dfl is given by 

r{cf>)=R{l + e^) (25) 

in polar coordinates {r,(j)). R is the radius of spiral at 
= and e is the deformation parameter determining 
relative size of the notch. The spiral-shaped billiard is 
fully chaotic, that is, there is no stable island at all due to 
its peculiar asymmetric property. We fix i? = 1 and con- 
sider two cases of deformation parameter, namely, weakly 
deformed case at e — 0.1 drawn in Fig. 5(a) and strongly 
deformed case at e = 0.3 depicted in Fig. 5(b). We take 
the spiral-shaped billiard as a model of non-convex geo- 
metric chaotic system and test the numerical procedure 
presented in Sec. II. 

The spiral-shaped microcavity laser has been firstly in- 
troduced by Chern et al. for obtaining unidirectional 
emission [42]. Afterwards, Lee et al. have found remark- 
able resonance patterns of the spiral-shaped dielectric mi- 
crocavity exhibiting strong localizations on a simple ge- 
ometric shape '43']. It looks like a clear counter-example 
of the conventional scar-theory in which the localized in- 
tensity patterns are appeared only on the corresponding 
classical unstable periodic orbits since the spiral has 
no simple-shaped periodic orbit, that is, all periodic or- 
bits must bounce the notch more than once [3]. Re- 
cently, Lee et al. have shown that above strongly local- 
ized resonance patterns can be approximated by linear 
combinations of nearly degenerated resonance modes of 
the circular cavity without any support from the classi- 
cal periodic orbits [i^ . Such recent research interests on 
the spiral-shaped microcavity also motivate the studies 
on spectral statistics of the spiral-shaped billiard. 

First we consider the weakly deformed spiral-shaped 
billiard at e = 0.1. We obtain 1500 accurate eigenvalue 

(a) (b) 




FIG. 5: The spiral-shaped billiards are shown on a same scale; 
(a) weakly deformed case for e = 0.1 and (b) strongly de- 
formed case for e = 0.3. 



spectra out of 21955 results calculated from the Ci con- 
tinuity quartic basis FEM. The (5„ fluctuates around zero 
in this range as shown in Fig. 6(a). In Figs. 6(b) and 
6(c), we present the results of I{s) and A3(L). We ex- 
pect that the spectral statistics are described by the GOE 
prediction due to the time reversal invariance as the case 
of cardioid billiard. However the results follow the GUE 
prediction rather than the GOE in Figs. 6(b) and 6(c). 

For the strongly deformed spiral-shaped billiard at 
e — 0.3, we attain the same 1500 accurate eigenvalues 
among 21934 calculated data (see the results of J„ in 
Fig. 7(a)). In contrast to the weakly deformed case, the 
results of I{s) and A3(L) are well described by the GOE 
expectation as one can show in Figs. 7(b) and 7(c). In 
result, different degrees of deformation causes quite dif- 
ferent spectral statistics and the unexpected GUE-like 
statistics are observed. 

In the literatures, there have been several reports that 
study the spectral statistics exhibiting unexpected GUE- 
like behavior in time reversal invariant systems, for exam- 
ples, the system with certain point symmetry (29l.l46j and 
the so-called Monza billiard possessing the property of 
unidirectional motion [20] . For these systems it has been 
known that the GUE-like spectral statistics have their 
origin in the degenerated eigenstates. However the GUE- 
like statistics of the weakly deformed spiral-shaped bil- 
liard at e = 0.1 cannot be explained by this reason, since 
the system has no degenerated eigenstates. Note that 
reasonable accounts for the unexpected spectral statis- 
tics of the weakly deformed spiral-shaped billiard are not 
feasible at present. But we would anticipate that all of 
above GUE-like spectral statistics can be understood in 
an unified principle. 



IV. CONCLUSION 

We present an efficient finite element method for cal- 
culating eigenvalues and eigenfunctions of quantum bil- 
liard systems. The Ci continuity quartic interpolation 
basis is considered. We show that the method provides 
accurate set of eigenvalues exceeding a thousand levels 
for any shape of quantum billiards on a personal com- 
puter. Comparison with the well-known Cq continuity 
quadratic basis FEM proves the efficiency of the applied 
method. The spectral statistics of the Robnik and the 
spiral-shaped billiards are studied and the unexpected 
GUE-like behaviors are observed. 

Note that we do not make use of the sparsity of ma- 
trices. We would expect that the generalized eigenvalue 
solving routine, which is optimized to sparse matrix, en- 
hances the efficiency of presented FEM. 
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FIG. 6: (Color online) For the weakly deformed spiral-shaped billiard ai e — 0.1; (a) The S„ up to 1500 level, (b) The cumulative 
level spacing distributions of the GOE, GUE, and the considered billiard are drawn by a black full, black dotted, and orange 
full line, respectively, (c) The spectral rigidities of the GOE, GUE, and the considered billiard are drawn by a black full, black 
dotted, and orange full line, respectively. 
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FIG. 7; (Color online) For the strongly deformed spiral-shaped billiard at e = 0.3; (a), (b), and (c) contain the same quantities 
of Fig. 6. 
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